function eq_mit = solve_eq_mit(sim, eq_SS, param, glob, options)

    diff = 10;
    %iter = 0;
    
    T_plot = 20;
    
    while (diff > 1e-7)
        
        sim_out  = solve_bf(sim, eq_SS, param, glob, options);
        
        diffC = max(max(abs(sim_out.C - sim.C)./abs(sim.C)));
        diffW = max(max(abs(sim_out.W - sim.W)./abs(sim.W)));
        diffSDF = max(max(abs(sim_out.sdf -sim.sdf)./abs(sim.sdf)));
        
        diff  = max([diffC diffW])
        
        Labor       =  sim_out.l;
        total_labor =  sim_out.L; sum(sim_out.mu.*Labor);
        
        figure(47)
        
        subplot(3,3,1)
        hold off
        plot([sim_out.C - 1]); drawnow; 
        hold on
        plot(sim.C-1); title('C'); drawnow;
  %      xlim([1 T_plot]);
        
        figure(47)

        subplot(3,3,2)
        hold off
        plot(sim_out.D - eq_SS.D); drawnow; title('D'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        
        subplot(3,3,3)
        hold off
        plot([total_labor]); title('Total labor demand'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        
        mass = sum(sim_out.mu);
        
        subplot(3,3,4)
        hold off
        plot(mass/mass(end)); title('Mass of firms'); drawnow;
        xlim([1 T_plot]);
        
        figure(47)
        subplot(3,3,5)
        hold off
        plot(sim_out.sdf); title('stochastic discount factor'); 
        hold on; plot(sim.sdf); drawnow;
%        xlim([1 T_plot]);

        figure(47)
        subplot(3,3,6)
        hold off
        plot(sim_out.ve); title('Value of entry'); drawnow;
        xlim([1 T_plot]);

        costs = sim_out.l;
        agg_prods = repmat(sim.A', glob.Nsf,1);
        prods = repmat(exp(glob.sf(:,2)), 1, options.T);
        prods = agg_prods.*prods;
        wages = reshape(sim_out.W,1,[]);
        wages = repmat(wages, glob.Nsf, 1);
        mkps  = sim_out.p./(wages./prods);
        
        tc    = sum(costs.*sim_out.mu);
        num    = sum(costs.*sim_out.mu.*mkps);
        
        mu_cw = num./tc;
        
        figure(47)
        subplot(3,3,7)
        hold off
        plot((mu_cw)); title('Cost-weighted markup'); drawnow;
        xlim([1 T_plot]);

        figure(47)
        subplot(3,3,8)
        hold off
        plot((sim_out.W)); 
        hold on
        plot((sim.W)); title('Wage');  drawnow;
 %       xlim([1 T_plot]);

%         % plot mass of entrants
         figure(47)
         subplot(3,3,9)
         hold off
         plot(sim_out.M);
         title('Mass of entrants');
         xlim([1 T_plot]);
         drawnow;
 
        sim.C =options.damp_mit.*sim.C + (1 - options.damp_mit).*sim_out.C;
        sim.W =options.damp_mit.*sim.W + (1 - options.damp_mit).*sim_out.W;
        
       sim.sdf = .99*sim.sdf + .01*sim_out.sdf;
%        sim.sdf = sim_out.sdf;           
%         if options.shock == 'tfp'
%             cEnew = param.ce*(sim_out.M./eq_SS.M).^(3/2);
%             sim.cE = .999*sim.cE + (1-.999)*cEnew;
%         end    
        
    end
    
    sim_out.entry = [];
    eq_mit = sim_out;
    
end